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Abstract 

The universality class of a phase transition is often determined by factors like dimensionality and 
inherent symmetry. We study the magnetic dipole system in which the ground-state symmetry and 
the underlying lattice structure are coupled to each other in an intricate way. A two-dimensional 
(2D) square-lattice system of magnetic dipoles undergoes an order-disorder phase transition be- 
longing to the 2D Ising universality class. According to Prakash and Henley [Phys. Rev. B 42, 
6572 (1990)], this can be related to the fourfold-symmetric ground states which suggests a similar- 
ity to the four-state clock model. Provided that this type of symmetry connection holds true, the 
magnetic dipoles on a honeycomb lattice, which possess sixfold-symmetric ground states, should 
exhibit a Kosterlitz-Thouless transition in accordance with the six-state clock model. This is ver- 
ified through numerical simulations in the present investigation. However, it is pointed out that 
this symmetry argument does not always apply, which suggests that factors other than symmetry 
can be decisive for the universality class of the magnetic dipole system. 

PACS numbers: 75.70.Ak,75.10.Hk,64.60.Cn 
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Understanding the physics in thin films is of practical importance since thin-film con- 
struction is used in manufacturing a variety of electronic and optical devices. One may 
note that magnetic thin films, in particular, play a key role in massive data storage applica- 
tions and that magnetic interactions in such films will become more prominent as magnetic 
moments are more densely integrated on such devices. A typical behavior of the magnetic 
property in thin films is switching between perpendicular and in-plane magnetization as 
the temperature T varies, which was first reported for Fe films.- On the other hand, for 
rare-earth compounds such as ErE^CusOg+z, the ordering of magnetic spins at T < 1.0K is 
essentially two-dimensional (2D).— For these materials, the exchange interaction is known to 
be relatively weak,- and the rare-earth ionic moments have been described as Ising or XY 
spins depending on anisotropy.- These compounds have also attracted attention as a suitable 
candidate to reveal the relation between superconductivity and magnetism. A model of a 
magnetic thin film is therefore a 2D lattice of magnets governed by the dipole interaction 
and confined to rotate on the plane of the lattice. Although it is straightforward to write 
down the corresponding Hamiltonian, understanding its physics is more cumbersome, for 
two reasons: the long-range character of the dipole interaction and its anisotropy. For a 
lattice of N planar magnets with the dipole interaction, the Hamiltonian is given as follows: 



where J(> 0) is a coupling constant and the summation runs over all the distinct spin pairs 
Si and Sj, residing at r» and rj, respectively. The distance between members of a spin 
pair is denoted as = \ri Since we employ the periodic-boundary condition, the 

displacement between r, and r 3 - is chosen as the one with the minimal distance among 
every possible pair of their periodic images. If more than two periodic images of a spin have 
the same minimal distance from another spin, we neglect the interaction between members 
of this spin pair to remove the ambiguity. As is clearly seen in Eq. (pQ), the interaction energy 
between Si and Sj decays as r^ 3 and it does not depend solely on their angular difference but 
also on their relative position, r^. In terms of numerical analysis, the long-range character 
imposes an 0(N 2 ) complexity within the simple Metropolis algorithm and the anisotropy 
puts an obstacle to developing an effective cluster algorithm. These problems have left the 
properties of the phase transition in this system largely inconclusive. On the one hand, 
this dipole lattice is related to the ice model,- which undergoes a phase transition keeping 
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its structural arrangement disordered.- The lack of a long-range order has been reported 
in the neutron-scattering experiments for some rare-earth compounds^ and also in artificial 
spin ice.- On the other hand, such disorder-preserving behavior apparently disagrees with 
theoretical and numerical predictions that the long-range order will be established in the 2D 
dipole lattice at low T.-~— Even if we accept the existence of an order-disorder transition 
at a critical temperature T c in the dipole system, numerical studies produce conflicting 
results: For the case of a square lattice, there has been reported a value of the critical 
exponent (3 = 0.19(4) for the staggered magnetization as well as 7 = 1.37(7) for the staggered 
susceptibility with the correlation- length exponent v = 1,— ^ where the numbers in the 
parentheses are numerical errors in the last digits. The staggered magnetization vector is 
defined as m = iV _1 Y2i ^ii where we specify the components of the vectors as = (xj, jji) 
and Si = (cos 8i, sin 6 i), respectively, and define gauge-transformed spins as 9 

^[(-l^cosM-lf'sinfli]. (2) 

We take the magnitude m — \m\ as a scalar magnetic-order parameter of this system. Us- 
ing the same observable, a recent study reported (3/v = 0.13(2) and v = 1.05(5) with the 
Metropolis algorithm.— This result indicates the 2D Ising universality class within errors, 
which is partially supported by another study on the Heisenberg dipole system.— How- 
ever, none of these match a renormalization-group calculation^ yielding an exponentially 
diverging correlation length £ with log£ ~ ~ T c and magnetization ~ £ -1 / 2 . To our 

knowledge, no theoretical explanation of the observed results has been satisfactorily pro- 
vided. 

In this brief report, we begin with the critical behavior of magnetic dipoles on the square 
lattice. The finite-size scaling result from L x L square lattices supports an order-disorder 
transition with the 2D Ising universality class, confirming the recent observation.— We 
then ask if this is related to the fact that the system possesses fourfold-symmetric ground 
states.- To test this further, we use the fact that the magnetic dipoles may have sixfold- 
symmetric ground states if put on the honeycomb lattice.- 6 This has been well established 
for the nearest-neighbor dipole interaction,— and we have numerically checked its validity 
for the long-range-interaction case as well. We find that the honeycomb-lattice lattice case 
exhibits a Kosterlitz-Thouless (KT) transition as implied by the analogy to the six-state 
clock model.— ^ 
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Our simulation strategy is as follows: we use the parallel tempering (PT) method that 
was devised to equilibrate glassy spin systems with very long relaxation times.— >^ Consider 
simulating two samples of a spin system in parallel at different inverse temperatures, 0i 
and /?2, respectively. If (3i < 02, the sample at /3i will explore a larger region in the phase 
space than will the sample at (3 2 . If the exploration happens to find a state with sufficiently 
low energy, we exchange Pi and 02 of the two samples so that the low-energy state can 
be pursued more deeply while the other sample begins a new exploration. Specifically, if 
energies of the two samples are denoted E\ and E 2 , respectively, the exchanging probability 
is given as P cx = min{l,exp [— (/3 2 — 0i)(E 1 — E 2 )]} to satisfy the detailed balance. One 
may easily extend this scheme to more than two samples, ranging over a broad temperature 
region, and run the samples simultaneously on parallel computing devices. We simulate 
each sample by the Metropolis algorithm, which means that the overall complexity is still 
0(N 2 ). Only its proportionality coefficient will be reduced by application of the PT method, 
but it is nevertheless a significant gain in practice especially when N is not too large. For 
L = 32 in Fig. GJa), for example, it usually takes 10 4 ~ 10 5 Monte Carlo steps for a simple 
Metropolis algorithm to equilibrate the system around T c from a random configuration, while 
it is enough to make a couple of exchange moves with running O(10 3 ) steps in between. We 
determine the difference in inverse temperature by observing overlaps of energy histograms. 
In the same example (L = 32), we have simultaneously simulated 31 inverse temperatures 
over [1.03, 1.63], each of which runs on an Intel Xeon quad-core L5420 CPU (2.5GHz). The 
CPU time spent for this size is about 6 x 10 2 h per temperature, meaning that the total 
CPU time to obtain the result for L = 32 roughly amounts to 2 x 10 h. 

Figure Ufa) shows Binder's cumulant U = 1 — | (m 4 ) / (m 2 ) 2 from the staggered magne- 
tization, where (■ • • ) means thermal average. Note that this cumulant is scaled to approach 
zero at high T and 1/2 at low T since (m 4 ) = 2 (m 2 ) 2 when the magnetization vector m 
has a 2D Gaussian distribution centered at the origin. From Fig. Ufa), the transition tem- 
perature is estimated as T c = 0.70(2) in units of J/ks, where ks is the Boltzmann constant. 
Since this is not an extremely precise estimation, it is hard to get critical exponents to a 
good precision. Instead, we may check consistency by assuming the 2D Ising universal- 
ity class. Let us plot the staggered magnetization and try a scaling collapse by the 2D 
Ising universality class, i.e., = 1/8 and v — 1 [Figs. Eta) andHfb)]. The best collapse 
is observed at T c = 0.71(1), which is in good agreement with the value of T c estimated 
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FIG. 1: (Color online) (a) Binder's cumulant of magnetic dipoles on L x L square lattices, (b) 
Cumulant ratio V^ 7-4 = (m^ / (m 2 ^ 2 obtained for the 2D square-lattice four-state clock model at 
the critical temperature T,? -4 = l/log(l + y/2). 
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FIG. 2: (Color online) Results on square lattices, (a) Staggered magnetization and its scaling 
collapse (b) according to the 2D Ising universality class, (c) Staggered susceptibility and its scaling 
collapse (d) in the same way. The 2D Ising critical exponents /3 = 1/8, v = 1, and 7 = 7/4 are 
used; and the critical temperature is estimated to be T c = 0.71(1) and 0.72(1) for (b) and (d), 
respectively. 

above via Binder's cumulant. Another important observable is the staggered susceptibility 
X = NT' 1 [(m 2 ) — (m) 2 ] plotted in Fig. El^c). Again the 2D Ising universality class with the 
corresponding critical exponent 7 = 7/4 is clearly consistent with the data provided that 
T c = 0.72(1) [Fig. EJ^d)]. One also notes that the susceptibility data imply a single transi- 
tion point where the critical fluctuation diverges in the thermodynamic limit, and that the 
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susceptibility remains finite below this divergence. This rules out the possibility of a KT 
transition. The consistent descriptions strongly suggest that the critical behavior indeed 
belongs to the 2D Ising universality class. 

It is obviously nontrivial that a continuous-spin system, complicated by the long-range 
character and anisotropy, nevertheless displays an Ising transition. The most plausible 
explanation is that the ground states of the dipole system on the square lattice possess a 
fourfold symmetry- since the same is true for the four-state clock model, which exhibits 
the 2D Ising criticality.— The value of the cumulant indeed tells us more than this simple 
symmetry argument: Figure Q^a) shows that V4 = (m 4 ) / (m 2 ) 2 = 2(1 — U) ~ 1.04(1). For 
the 2D Ising model, on the other hand, the value of this quantity is given as V^ Ismg = 1 + e, 
where e = 0.167923(5).™ Since the four-state clock model is equivalent to two independent 
Ising systems A and B with the temperature rescaled,— we can denote it as A <g> B and 
consider its magnetic-order parameter m = \J{m 2 A + m 2 3 )/2 where and m# correspond 
to magnetizations of the independent Ising systems. This leads to the cumulant value of the 
four-state clock model as follows: 



which agrees well with our Monte Carlo calculation V£~ 4 = 1.0840(2) at T«= 4 = l/log(l + 
V%) [Fig. 00(b)]. By the same analogy, the cumulant value for the magnetic dipoles sug- 
gests that V4 = l + e/4 w 1.042, which could be observed in a combination of four 
independent Ising systems denoted as A, B, C, and D, respectively, or two independent 
four-state clock models A <g> B and C <8> D, where the total magnetization is equivalent to 
m = (m\ + m 2 B + m? c + m 2 D ) / 4. 

In order to test the symmetry argument further, we investigate the honeycomb lattice, 
since the nearest-neighbor dipole interaction in this case leads to sixfold ground states.— In 
addition, since it is known that the g-state clock model with the cosine interaction undergoes 
a KT transition when q > 6,—— the symmetry argument suggests that the transition for a 
honeycomb lattice should be of a KT type. 

The honeycomb lattice that we use in this work is shown in Fig. [3](a) where the system size 
is iV = 2L 2 for a given length scale L. We choose L as multiples of 6 by taking the ground- 
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FIG. 3: (a) Honeycomb lattice of size L = 6 with the periodic-boundary condition. Coordinates 
are presented as (u, w) for the gauge transformation in Table HI (b) Distribution of staggered 
magnetization m = (m x ,my), showing the sixfold symmetry, taken at T = 0.1 with L = 12. 

state configurations^ into consideration. An appropriate gauge transformation of the spin 
angle 6 at site (u, w) in the same spirit of Eq. ([2]) is tabulated in Table HI The staggered 
magnetization is defined as m = iV _1 ^(cos#', sin^') with magnitude m = \m\. As in 
the square-lattice case, each chosen spin interacts with a set of other spins that have well- 
defined minimal distances to the chosen spin. We furthermore require this interacting set to 
have the inherent symmetry of the honeycomb lattice; i.e., the interacting set should be left 
invariant under the rotation by ±27r/3 around the chosen spin so that the six ground states 
are equally probable in the system [Fig. |3](b)] . The numerical results are again obtained by 
the PT method and plotted in Figs.|4](a)-|4](d), where Binder's cumulant U and the staggered 
susceptibility x are given in the same way as above. In order to examine the low-T phase, 
we have the inverse temperature range over 0(1) to O(10) while keeping the overlaps in 
the energy histograms. For L = 24, for example, we ran 127 inverse temperatures over 
[1.5, 10.3] in parallel, spending about 10 2 CPU hours per each. Around T 0.5, one finds a 
size merging of U, together with a divergence in x, which are characteristic signatures of the 
KT transition. These observations imply the existence of a quasicritical phase below T w 0.5 
where the correlation length diverges. The KT picture also predicts a scaling collapse of x vs 
U, such that x ~ L 2 ~ v f(U) where / is a certain scaling function and r\ = 1/4.— This method 
provides a piece of information about the universality class even without precise knowledge 
of the transition temperature. In spite of the small system sizes, the KT scaling exponent 
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TABLE I: Gauge transformation for the honeycomb lattice, where (u, w) represents the coordinates 
as shown in Fig. [3ta). 
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77 = 1/4 nevertheless gives a scaling collapse consistent with a KT transition as shown in 
the inset of Fig. Etc). We also note an additional tiny yet systematic size dependence of U 
below T = 0.2 shown in the inset of Fig. HJ^b), which possibly indicates that the staggered 
magnetization freezes into the sixfold symmetry [compare Fig. [3(b)] . Unlike in the six- 
state clock model, however, this freezing is not accompanied by any peak in specific heat 
[Fig.HJd)], which implies that the quasicritical phase is not isotropic either but should reflect 
the sixfold symmetry at least in part. It is currently under investigation whether the U(l) 
symmetry in the disordered phase gets broken exactly at the same temperature where the 
KT transition occurs. 

In summary, we have confirmed that the XY-type magnetic dipoles on the square lattice 
exhibit the 2D Ising criticality and that this can be related to a symmetry similarity in 
the four-state clock model. This symmetry connection is further supported by the study 
of the honeycomb lattice, where the ground states have sixfold symmetry; and the system 
behaves similarly to the six-state clock model exhibiting a KT transition. However, the 
symmetry argument is not always the decisive factor: the transition is KT-like for the square 
lattice with nearest-neighbor interaction in spite of the fact that the symmetry remains the 
same. 24 One may also note that the long-range order at low T is also absent in experiments 
with squarelike structures.- However, the factors that supersede the symmetry argument in 
deciding the critical universality class remain to be elucidated. 
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FIG. 4: (Color online) Results on honeycomb lattices, (a) Magnitude of the staggered magne- 
tization and (b) Binder's cumulant, with a closer view of the low-T part (inset), (c) Staggered 
susceptibility \ an d the scaling collapse of \ ~ L 2 ~ v f(U) with 77 = 1/4. (d) Specific heat does not 
show any clear double-peak structure. 
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